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1.  INTRODUCTION 


Eigenanalysis  of  temperature  and  wind  profiles  often  reveals  distinct  resonances  at 
mesoscale  horizontal  wavelengths  which  are  on  the  order  of  2  km  to  30  km.  An  example 
of  such  an  analysis  is  described  in  Appendix  A.  These  analyses  suggest  using  linear 
methods  that  can  identify  and  depict  mesoscale  resonant  modes  such  as  mountain  waves. 
Fourier  transform  analysis  methods,  when  properly  posed,  appear  well-suited  for  this 
purpose. 

The  strengths  of  Fourier  models  apply  to  several  areas.  Gravity  wave  propagation 
methods  using  Fourier  analysis  can  run  faster  on  computers  than  mesoscale  numerical 
prediction  (NWP)  models  with  comparable  spatial  resolution.  This  allows  use  of  small 
computers  in  many  applications.  The  linear  assumptions  inherent  in  Fourier  analysis 
require  only  a  horizontally-uniform  background  state  which  can  be  provided  from 
relatively  coarse  numerical  model  output  or  rawinsonde  data.  The  quasi-analytical 
solutions  avoid  numerical  methods  with  their  attendant  error.  Fourier  methods  use 
complex  arithmetic  which  allows  both  modulus  and  phase  information  in  the  wave 
function.  These  methods  also  permit  continuous  functional  behavior  in  the  vertical. 

Fourier  methods  used  for  mountain  wave  analysis  also  bear  some  resemblance  to  wave 
propagation  methods  employed  in  the  physical  sciences  and  engineering  for  determining 
wave  dispersion  and  attenuation.  This  includes  fields  as  distinct  as  oceanography, 
seismology,  optics,  and  radio  science. 

There  are  notable  limitations  to  Fourier  analysis  methods.  A  confined  horizontal  domain 
must  typically  be  considered  for  analysis  given  the  assumption  of  horizontally  uniform 
background  flow.  Also,  Fourier  models  cannot  treat  non-linear  behavior,  such  as 
breaking  waves,  where  small-amplitude  wave  perturbations  interact  with  the  larger  scale 
uniform  basic  state.  This  inability  to  handle  non-linear  behavior  may  be  its  greatest 
limitation. 


2.  COMMON  TRAITS  OF  FOURIER  MOUNTAIN  WAVE  MODELS 


Common  traits  of  Fourier  mountain  wave  models  include  the  use  of  the  Boussinesq 
approximation  where  the  density  is  constant  except  for  terms  involving  the  external  force 
(e.g.  gravity)  in  the  equations  of  motion  (Chandrasekhar,  1961).  The  Boussinesq 
approximation  allows  no  acoustic  waves  and  pennits  only  shallow  oscillations  relative  to 
the  atmospheric  scale  height.  The  anelastic  approximation,  which  allows  for  deep 
motion,  is  sometimes  used  instead  and  typically  results  in  the  amplitude  being  scaled  by 
the  inverse  square  root  of  the  density. 


The  Fourier  models  represent  solutions  to  elliptic  partial  differential  equations,  expressed 
in  terms  of  height  or  vertical  velocity  perturbation.  These  equations  are  derived  from  the 
governing  system  of  three-dimensional  (3-D)  linearized  partial  differential  equations  for 
conservation  of  mass  and  momentum,  and  written  in  the  complex  wave  number  domain 
in  the  form: 

i^(fj2 it* +  ^k*  +  ky2)(N2  -a2)v  =  0  ( m  s~2  rad) 

Here  77  is  the  three-dimensional  (3-D)  complex  height  perturbation  field,  kx  and  kv  are  the 

horizontal  wave  numbers,  N  is  the  Brunt-Vaisala  frequency,  z  is  the  vertical  coordinate, 
and  er  is  the  intrinsic  frequency.  Coriolis  acceleration  is  typically  neglected  for  the 
characteristic  time  scales  involved  with  mountain  waves  (~  30  minutes)  and  a  steady- 
state  background  atmosphere  is  assumed  in  accord  with  linear  theory.  The  models  are 
“dry”  in  that  no  moisture  effects  including  latent  heat  release  are  accounted  for  in  the 
equations.  Very  detailed  terrain  can  be  included  in  the  two-dimensional  lower  boundary 
with  its  detail  limited  only  by  the  horizontal  grid  spacing.  However,  in  accord  with  the 
linear  assumption,  the  maximum  value  for  77  at  the  lower  boundary  must  be  much  less 
than  the  total  model  depth. 

The  specification  of  the  three-dimensional  vertical  wave  number  or  phase  function 
relation  in  the  mountain  wave  model  is  critical.  Vertically-propagating  locally-plane 
waves  are  assumed.  For  vertically  propagating  wave  modes,  a  general  solution  of  the 
elliptic  equations  can  be  posed  as: 

n(kx,ky,z)  =  h  (kx ,  ky  )-g(kx,ky,z )  exp  [#  (kx  z)]  (m3/rad)  ( 1 ) 

where  the  complex  phase  function^  is  defined  in  integral  form  as 

< i(kx,ky,z)  =  -^m(kx,ky,z)dz  (rad)  (2) 

Here  m  is  the  vertical  wave  number  or  dispersion  relation,  h  is  the  terrain  lower 
boundary,  77  is  the  wave  height  perturbation  at  some  altitude  z,  and  g  is  a  density-  or 
altitude-dependent  function  related  to  the  vertical  wave  flux.  The  complex  variables  used 
in  this  expression  are  formulated  in  the  wave  number  or  transform  domain  as  represented 
by  the  hatted  symbol  A ,  where  i  =  4-i  .  Specific  forms  of  the  dispersion  relation  m  and 
function  g  for  each  model  will  be  discussed  in  more  detail  below. 

Some  general  relations  can  be  quantified  for  the  inclusion  of  phase  and  amplitude 
modulation  in  the  Fourier  models  relative  to  wave  dispersion.  The  vertical  wave  number 
777  or  dispersion  relation  is  expressed  as 

777  =  ^[k2x  +  kt2)(V2  -cr2)/crj^  (rad  m'1)  (3) 

where  N  is  the  Brunt-Vaisala  frequency  and  a  is  the  intrinsic  frequency  of  a  fluid  parcel 
moving  through  a  stationary  wave  field.  The  Brunt-Vaisala  frequency  is  the  square  root 
of  the  product  of  the  vertical  potential  temperature  gradient  dO /dz  (static  stability)  and 
the  ratio  of  acceleration  due  to  gravity  and  potential  temperature.  Specifically, 


(4) 


N  = 


g  dd 
6  dz 


(rad  s'1) 


The  Brunt-Vaisala  frequency  N  or  its  square  N2  is  also  referred  to  as  “static  stability”. 
Here,  the  over  bar  denotes  the  basic  state  or  horizontally-uniform  background  value. 
For  two-dimensional  horizontal  wind  flow,  the  intrinsic  frequency  is  defined  as 


a  =  -ky-kV 


(s'1) 


(5) 


where  U  and  V  are  the  west-to-east  and  south-to-north  background  wind  speeds, 
respectively.  For  the  dispersion  relation  used  throughout  this  report,  the  wind  curvature 
1 

terms  such  as  — - are  ignored  since  U  and  V  are  considered  to  be  slowly  varying  with 

U  dz 

height. 


Certain  relations  between  the  wave  dispersion  and  vertical  wave  propagation 
characteristics  exist  given  homogeneous  background  atmospheric  conditions  over  a 

rz 

propagation  layer  where  the  phase  angle  </>  =  -\  m  dz  for  a  vertically  propagating  wave.  If 

N2  -  a2  >  0 ,  then  m  is  real,  the  solution  is  trigonometric,  and  the  phase  is  modulated.  If 
N2  -  a2  »  0 ,  the  wave  solution  is  non-dispersive  and  vertically  propagating  and  m  is 
nearly  independent  of  cr  in  the  two-dimensional  limit.  If  N2  -  a2 ,  this  is  referred  to  as 
the  turning  point  height  or  buoyancy  frequency  turning  point  which  will  result  in  the 
trapping  of  waves.  The  role  of  the  turning  point  height  will  be  discussed  in  more  detail 
later  in  this  report.  If  N2  -  a2  <0 ,  then  m  is  imaginary,  the  solution  is  exponential,  and 
reflection  of  wave  energy  occurs  below  the  level.  If  a2  -  0 ,  a  singularity  exists  in  Eq.  (3) 
which  will  result  in  a  “critical  layer”  where  wave  absorption  will  occur  with  strong 
attenuation  above  the  level  of  analysis.  (Smith,  2002) 


3.  EVALUATION  OF  THE  TWO-DIMENSIONAL  FOURIER 
TRANSFORM 

The  mountain  wave  model  equations  at  a  given  height  are  solved  in  the  horizontal  two- 
dimensional  Fourier  transform  domain  for  a  function  /  ( x,y,z ) : 

F(kx,ky,z)  =  J J(x,y,z)  exp[z'(krx  +  kyy)\dxdy  (6) 

using  the  Cartesian  coordinate  system  (x,  y,  z)  where  kx  and  k  are  horizontal  wave 

numbers  in  x  (west-to-east)  and y  (south-to-north)  directions,  respectively.  Working  in 
the  complex  transform  domain  allows  for  both  modulus  and  phase  information  in  the 
solution.  These  solutions  are  recovered  for  the  spatial  (and  real)  domain  using  the 
inverse  Fourier  transform. 

In  the  numerical  evaluation  of  the  double  Fourier  transform  such  as  Eq.  (6),  it  is  most 
efficient  to  express  this  function  as  two  successive  one-dimensional  transforms. 

Following  Sneddon  (1951), 


F(kx ,k)  =  -^-\  j  f(x, y ) e‘k-xe'k'ydx dy 
=  i —  f  i —  f  f(x,y)e,k'xdx  e'k,ydy 

yin  J“co|_V2/z' J_c0 

=  -7=  j ”  7 (7 , y)ek,ydy  =  F(kx , k  ) 
din  J-co 

This  results  in  2  n  operations  per  transform  instead  of  n4  operations  when  the  integrand 
of  Eq.  (6)  is  treated  as  a  four-dimensional  function.  Here  n  represents  the  number  of 
functional  data  points  to  be  transformed. 

Further  efficiency  is  obtained  by  using  the  fast  Fourier  transform  (FFT)  in  lieu  of  the 
mathematical  transform  shown  above.  The  fast  Fourier  transform  is  an  algorithm 
specially  designed  for  efficient  machine  calculation.  The  number  of  FFT  operations  for 
the  one  dimensional  (1-D)  Fourier  transform  is  proportional  to  n  log  n  rather  than  n  for 
the  direct  mathematical  form.  The  one  dimensional  FFT  algorithm  listed  in  Brigham 
(1974)  is  used  with  the  successive  1-D  transform  approach. 

4.  PRINCIPAL  FOURIER  MOUNTAIN  WAVE  MODELS 

Two  model  approaches  for  Fourier  mountain  wave  analysis  have  been  considered.  These 
include  the  Smith  (Yale  University)  Fourier  Three  Layer  model  and  the  Broutman  (Naval 
Research  Laboratory)  Mountain  Wave  Forecast  Model. 

The  Smith- Yale  Three-Layer  model  is  based  on  a  theoretical  foundation  that  has  been 
carefully  developed  since  1980.  As  its  name  implies,  it  uses  a  three  layer  profile  of 
constant  background  velocity  and  static  stability  within  a  layer  and,  in  general,  does  not 
pennit  wind  turning  with  height.  The  model  has  been  validated  during  field  tests 
including  the  Mesoscale  Alpine  Program  of  1999  (Smith  et  al,  2002). 

The  Broutman  model,  which  is  also  known  as  the  Mountain  Wave  Forecast  Model 
version  3  (MWFM-3),  uses  continuous  vertical  profiles  of  wind  shear  and  static  stability. 
A  transient  or  time-dependent  solution  is  also  available  for  trapped  waves. 

Both  Fourier  models  were  checked  by  their  authors  against  numerical  models  integrated 
in  time  to  a  steady-state  solution  using  the  same  basic  state  and  lower  boundary.  The 
Smith-Yale  model  was  checked  against  the  COAMPS  numerical  model  of  the  Naval 
Research  Laboratory  (NRL)  and  the  MWFM-3  model  was  checked  against  the  Lipps  and 
Hemler  (1982)  model  and  the  Weather  Research  and  Forecasting  (WRF)  model.  Some 
Fourier  and  numerical  model  comparisons  are  described  briefly  later  in  this  report. 

5.  SMITH  FOURIER  THREE-LAYER  MODEL 

A  two-dimensional  Fourier  transform  for  the  general  solution  of  the  elliptic  partial 
differential  equation  for  vertical  displacement  77  (in  meters)  is  (following  Smith,  2002) 


(m3  rad'1) 


(7) 


fj(k,l,z )  =  A1  exp(fmiz)  +  Bi  exp z) 
l  =  J  77 (x, j,z)exp^/^vx  +  kyy^\dxdy 

where  the  integers  i=  1,  2,  3  denote  three  layers  of  constant  wind  velocity  and  stability, 
Ai(kx,k\  and  Bj(kx,ky )  are  complex  coefficients  defined  below,  mi[kx,k v)  is  the 

complex  dispersion  relation  or  vertical  wavenumber,  and  z  is  the  height  (in  meters)  above 
the  model  lower  boundary.  The  vertical  wavenumber  m,  in  a  layer  i  for  non-hydrostatic 
flow  is  expressed  as 

m:  =  [k]  +  k;  )(iV,2  -  cx,2  )/cr,2  (rad  m'2)  (8) 
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where  kx  =  — ,  kr  =  —  (rad  /  m)  are  the  horizontal  wave  numbers  associated  with  the 

horizontal  wavelengths  Lx  and  Ly  respectively,  Nt  is  the  layer  Brunt- Vaisala  frequency 
(rad  sec'1),  and  cr.  =  kxUi  +  k  Vt  is  the  layer  intrinsic  frequency.  The  intrinsic  frequency  cr 

is  the  frequency  felt  by  a  parcel  of  fluid  moving  through  the  stationary  wave  field.  The 
vertical  wave  number  and  Fourier  coefficients  associated  with  a  particular  layer  i  will  be 
selected  according  to  the  height  or  layer  where  z  is  defined. 

The  Fourier  amplitude  coefficients  T.  (kx,ky  J ,  B  ( kx ,  ky  j  are  computed  using  the  transform 
of  the  terrain  i(kx,k\{  m3  radian'1)  and  the  dispersive  parameters  that  define  each  layer 
(/rtj.jtr^z,.) .  The  expressions  for  three  layers  in  the  Smith  model  are: 

A3=h/(FAE  +  q-FBE) 

A ,  =  FAE  ■  A3 
B  .t  =  FBE  ■  4 
4  -  FCE  ■  4 
B  -,  -  FDE  ■  4 

where 

FCE  =  (1/2) [l  +  R32  ]  exp (im3z2  -  im2z2 ) 

FDE  =  (l/2)[l-42]exP(”M3z2  +  im2z2 ) 

FAE  =  (l/4)  [(l  +  R2l )  (l  +  R32 )  J  exp  ( i(m2  -  mx  )zj  +  i  ( m3  -  m2  )z2 ) 

+(l/4)[(l-4i)(l-42)]exp(-d>2  +  m\)z\  +i(m  3  +  w2)z2) 

FBE  =  (1/4)  [(1  -  R2l  )(1  +  R32  )]exp  (i(m2  +  mx )z,  +  i(m3  -  mn  )z2 ) 

+(1/4)[(1  +  4i)(l-42)]exp(-7(m2  -  mx  )zj  +  i(m3  +  »i2  )z2 ) 


and  where 


^32  =m3aj/m2  u\ 

Rn  =  m2c?llml  of 

Zj  here  is  defined  as  the  top  of  layer  1  and  z2  is  the  top  of  layer  2.  All  quantities  are 
treated  as  complex.  S3  is  set  equal  to  zero  for  the  top  layer  using  the  radiation  condition 
at  the  upper  boundary. 

To  give  the  model  flexibility,  a  reflection  coefficient  q  ,  different  from  unity,  is 
introduced  for  the  lower  boundary  to  represent  partial  absorption  of  down-going  waves. 
Dissipation  of  the  wave  by  boundary  layer  turbulence  or  by  critical  layer  absorption  at  the 
lower  boundary  can  be  parameterized  by  setting  q  <  1 .  When  atmospheric  conditions  are 
uniform  with  height,  then  R2l  =  Z?32  =  1  and  the  down-going  wave  amplitudes 
are  5,  =  B2  =  0  .  The  base  level  in  the  calculation  is  taken  to  be  the  reference  height  zref  or 

the  top  of  the  pre-defined  stagnant  layer.  The  linearized  lower  boundary  condition  is  then 
the  height  h  of  the  terrain  above  the  reference  height  zref : 

ri{x,y,ti)  =  h{x,y)-zref  (9) 

Equation  (7)  is  solved  for  fj  for  any  level  z  using  the  Fourier  amplitude  coefficients  and 
dispersion  relations  that  correspond  to  the  layer  where  z  is  located.  If  z  resides  on  a  layer 
boundary,  the  coefficients  in  Eq.  (7)  are  assigned  for  the  layer  below  z .  Diagnostic 
expressions  for  the  perturbation  vertical  velocity  amplitudes  w  and  u  are  calculated  in  the 
Fourier  domain  using  77 ,  the  intrinsic  frequency,  and  certain  wave  number  relations. 
These  velocity  tenns  are  expressed  as: 

w-iaq  (m3  s"1)  (10) 

„  kNw  ,  3  _k 

u  = - ~p~ - -  (m  s  )  (11) 

*klf{k,u+k,y) 


An  inverse  Fourier  transform  is  then  perfonned  to  display  these  variables  in  the  real 
domain,  in  the  form: 

i j(x,y,z)  -  ^  J  J)[kx,ky,z)e~lKx e~lkyy dkxdky  (m) 

The  perturbation  temperature  O'  at  a  given  point  can  be  calculated  using  77  and  the  vertical 


SO 

gradient  of  mean  potential  temperature  — : 

dz 


„  dO 

0  =77  — 

dz 


(deg  K)  (12) 


5.1  Computational  Aspects  of  the  Smith  Three-Layer  Model 

Solutions  at  a  given  altitude  are  obtained  for  a  1024  x  1024  domain  of  horizontal  grid 
points.  This  corresponds  to  a  2P  dimension  using  the  base  2  FFT  for  p  =  10.  The 


horizontal  grid  point  spacing  Ax,  Ay  is  1  km  resulting  in  a  bandwidth  or  wave  number 
differential  Akx  =  A ky  =  2;r/(«Ax)  =  6.136  x  10'6  rad/m. 

For  the  terrain  conditions,  30  arc-second  terrain  data  from  the  NOAA  GLOBE  database 
(Hastings  and  Dunbar,  1998)  are  used  to  construct  the  lower  terrain  boundary 
height  h(x,y) .  The  GLOBE  database  terrain  heights  are  defined  as  a  function  of  latitude 

and  longitude,  which  corresponds  to  less  than  one  km  grid  spacing  for  both  latitudinal 
and  longitudinal  directions  when  it  is  mapped  to  a  Cartesian  grid.  The  latitudinal  spacing 
is  more  or  less  constant  with  latitude  while  the  longitudinal  grid  spacing  decreases  with 
increasing  poleward  latitude.  For  the  model’s  purpose,  the  GLOBE  terrain  height  data 
are  interpolated  to  one  km  increments  for  the  horizontal  grid. 

Equation  (7)  is  evaluated  at  any  wave  number  grid  point  at  a  given  level  except  for  those 
where  critical  levels  or  singularities  occur  in  the  dispersion  relation  (i.e.,  a  -  0).  These 
grid  points  are  skipped  and  are  not  included  in  the  transform  solution.  Also,  both  a 
reference  level  and  lower  boundary  reflection  coefficient  q  must  be  defined.  The 
reference  level  serves  to  reduce  the  lower  model  boundary  terrain  height  by  the  height  of 
the  reference  level  above  sea  level.  This  reference  level  typically  corresponds  to  the 
height  of  the  boundary  or  “stagnant”  layer.  The  lower  boundary  reflection  coefficient 
appears  in  the  expression  for  the  3ld  layer  coefficient  i  =  3  and  is  set  to  q  <  1 . 


5.2  Practical  Considerations 

The  terrain  area  h(x,y)  for  the  lower  boundary  condition  is  first  padded  by  “zeroes” 

around  its  horizontal  boundary  for  dimensional  expansion  of  the  computation  domain 
necessary  to  meet  numerical  Fourier  transform  requirements.  Typically,  the  expansion 
satisfies  2P  grid  points  in  a  single  dimension  for  use  in  the  FFT  algorithm  where  p  is  an 
integer.  The  zero-padding,  by  definition,  also  minimizes  truncation  error  of  the  discrete 
Fourier  transform.  However,  truncation  leads  to  a  “wrap-around”  effect,  associated  with 
periodic  boundary  conditions,  where  solutions  occurring  outside  of  the  terrain  area  on 
one  side  of  the  computational  domain  will  re-appear  on  the  opposite  side  of  the  domain. 
For  certain  classes  of  terrain,  an  exponential  decay  with  distance  of  terrain  height  from 
the  high  elevation  point  of  the  grid  is  introduced.  This  can  smooth  the  joining  of  terrain 
boundaries  and  padding  to  minimize  spurious  waves  introduced  by  discontinuities 
between  the  terrain  and  padded  areas. 

These  terrain  padding  and  transfonn  operation  procedures  are  illustrated  in  Figure  1. 
After  subtracting  the  reference  level  from  the  terrain  height,  the  terrain  field  is  padded 
with  zeroes  to  minimize  “wrap-around”  associated  with  periodic  function  truncation.  If 
appropriate,  an  exponential  decay  of  terrain  height  from  the  location  of  maximum  height 
location  may  be  introduced.  This  yields  the  field  of  h[x,y )  depicted  at  the  top  of  Figure 
1.  The  two-dimensional  FFT  of  the  h(x,y)  field  yields  the  complex  transform  field 
H(kx,  ky )  =  ^  /7r  ( kx ,  ky ) ,  /7,(  kx ,  ky )  J  ,  with  Hr  and  Hi  representing  the  real  and  imaginary 


transform  components  respectively,  which  is  depicted  in  the  middle  figure.  This  figure 
represents  the  operational  wave  number  domain  where  the  model  solutions  are  obtained. 
The  inverse  transfonn  of H(kx,ky)  yields  (approximately)  the  original  function  h  (x,y) 

which  is  shown  at  the  bottom  of  Figure  1. 


h(x,y) 


Figure  1.  The  Fourier  transform  pair  for  terrain  height  in  the  real  domain  h(x,y'j  and  the 
complex  terrain  H(kx,ky )  of  the  wave  number  domain. 


5.3  Examples  of  Smith  Three-Layer  Model  Results  using  Theoretical 
and  Real-World  Terrain  and  Basic  State 

5.3.1  THEORETICAL  RESULTS 

To  check  the  validity  of  the  Smith  Three  Layer  model  implementation,  four  theoretical 
cases  described  in  Smith’s  (2002)  monograph  were  evaluated  and  compared  with  his 
published  results.  The  calculations  matched  the  published  results  and  also  serve  to 
illustrate  wave  behavior  for  different  combinations  of  stability,  vertical  shear,  and  terrain 
shape  and  terrain  orientation.  These  four  cases  included  (1)  an  axisymmetric  hill  with 
uniform  wind  velocity  and  stability,  (2)  a  ridge  line  perpendicular  to  the  mean  flow 
direction  with  uniform  wind  velocity  and  stability,  (3)  an  axisymmetric  hill  and  ridge  line 
with  decreasing  Scorer  parameter  (N/U)  with  height,  and  (4)  the  ridge  line  oriented  45°  to 
the  mean  flow  with  decreasing  Scorer  parameter  with  height.  The  wind  direction  is 
constant  with  altitude  for  all  cases.  The  inputs  for  the  theoretical  calculations  are  listed  in 
Table  1. 


For  the  hill  shape,  the  equation  representing  the  terrain  height  is 

z(x,y)  =  A,ax  exp|~-(x/a)2  - (y/b)2 ]  (m) 


where  a  and  b  are,  respectively,  the  minor  and  major  axes  of  the  elliptical  mountain 
shape.  For  the  axisymmetric  hill  shape,  a  =  b.  In  Figures  2  and  3,  the  planform  analysis 
height  is  5  km,  while  the  planform  analysis  level  is  at  2  km  in  Figures  4  and  5.  The 
vertical  cross-section  is  valid  for  the  red  line  positioned  at  y  =  0.  The  results  are 
displayed  in  terms  of  height  deviation  A z  from  the  mean  or  background  height  with  /*max  = 
1000  meters. 


Table  1.  Model  Inputs  for  Theoretical  Calculations  Using  Three  Layer  Model 


Terrain  ( a ,  b ) 
hill/ridse 

Stability 

N(  s'1) _ 

Wind  Speed 

U  (m/s) 

Scorer  Parameter 
N/U  Cm  '  ) 

Absorption 

Fig.  2 

10  km 

0.01 

10 

0.001 

none 

Fig.  3 

10/50 

0.01 

10 

0.001 

none 

Fig.  4a 

2  km 

0.012/0.008 

10/22 

0.0012/0.00036 

q  =  0.93 

Fig.  4b 

2/50 

0.012/0.008 

10/22 

0.0012/0.00036 

q  =  0.90 

Fig.  5a 

2/50 

0.012/0.008 

8/22 

0.0015/0.00036 

q  =  0.93 

Fig.  5b 

2/50/45° 

0.012/0.008 

10/22 

0.0012/0.00036 

q  =  0.90 

Figures  2  and  3  depict  examples  of  terrain-induced  waves  in  a  hydrostatic  environment 
since  little  energy  will  be  put  into  wavelengths  less  than  the  hill  radius.  In  Figure  2,  the 
energy  is  confined  to  a  downstream  parabolic  region  while  for  the  ridge  case  in  Figure  3, 
the  energy  is  found  only  in  the  vicinity  of  the  ridge.  In  both  cases,  the  wave  phase  tilts 
upstream  into  the  wind.  In  Figure  4,  the  waves  generated  by  the  axisymmetric  hill  and 
ridge  line  are  trapped  by  the  change  in  wind  speed  and  stability  at  2  km,  resulting  in 
transverse  wave  modes  (wave  bands  perpendicular  to  the  shear  vector)  as  well  as  a 
diverging  parabolic  wave  for  the  axisymmetric  hill.  In  Figure  5,  the  wind  speed  and 


stability  changes  are  introduced  at  a  higher  altitude  associated  with  a  deeper  stable  layer. 
This  results  in  a  short  wavelength  fundamental  mode  and  a  second,  longer  wavelength 
mode.  In  Figure  5b,  the  ridge  is  rotated  clockwise  45°  which  results  in  two  modes  where  the 
longer  wavelength  mode  is  rotated  farther  to  the  north  than  the  shorter  wavelength  mode. 


axisymmetric  hill  (10  km) 


Az  at  5  km 


-100  -50 


50  200 


-  uniform  profile 

-  hydrostatic 

-  constant  wind  speed  and  stability 

-  total  absorption  at  lower  boundary 


s  (km) 

Figure  2.  Gaussian  hill  profile  (Witch-of-Agnesi)  with  constant  wind  speed  and  stability.  Total 
absorption  at  lower  boundary. 


ridge  (10  x  50  km) 


Az  at  5  km 


'-100  -50  0  50  1  00  1  50  200 


-  uniform  profile 

-  hydrostatic 

-  constant  wind  speed  and  stability 

-  total  absorption  at  lower  boundary 


Figure  3.  Ridge  profile  with  constant  wind  speed  and  stability.  Total  absorption  at  lower 
boundary. 


axisymmetric  hill  (2  km) 


ridge  (2  km  x  50  km) 


50- 


Az  at  2  km 


-100  -50  0  50  100  1  50  200 


-  two-layer  Scorer  number  (N/U)  profile 

-  N/U  decreasing  with  height  (2.4  km) 

-  evanescent,  trapped  waves 


Figure  4.  Two-layer  wind  speed  and  stability  with  reflection  at  lower  boundary, 
a.  axisymmetric  hill ,  q  =0.93  (reflectivity  factor)  b.  ridge  (2  x  50  km),  q=0.90 


ridge  (2  km  x  50  km) 


100- 


-100-^ 
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Az  at  2  km 


ridge  (2  km  x  50  km)  rotated  by  45° 


-  two-layer  Scorer  number  profile 

-  lower  layer  depth  increased  to  4  km 

-  2  transverse  modes  (  Ax  ~  5 , 14  km) 


Figure  5.  Two-layer  wind  speed  and  stability  with  reflection  at  lower  boundary. 

a.  ridge  (2  x  50  km) ,  q  =0.93  (reflectivity  factor)  b.  ridge  (2  x  50  km)  rotated  by  45°,  q=0.90 


5.3.2  REAL-WORLD  TERRAIN  AND  BASIC  STATE 


Calculations  using  the  Smith  Mountain  Three  Layer  Wave  Model  are  illustrated  using 
real-world  terrain  over  the  French  Alps  for  a  case  study  described  by  Smith  et  al  (2002) 
during  the  1999  Mesoscale  Alpine  Programme  (MAP)  campaign.  This  is  the  same  terrain 
field  depicted  in  the  terrain  padding  example  described  earlier.  The  terrain  field  is 
centered  near  latitude  46°  N,  longitude  7°  E.  A  three  layer  model  profile  of  wind  velocity 
and  static  stability  is  used  which  is  based  on  rawinsonde  and  dropsonde  measurements 
taken  at  1200  UTC  on  2  November  1999.  Given  a  stagnant  layer  top  of  2.5  km  above 
mean  sea  level  (MSL)  and  a  reflectivity  factor  q  =  0.90,  Table  2  lists  the  layers  and  their 
mean  static  stability  and  mean  wind  velocity  characteristics  used  for  the  model  input. 


Table  2.  Three-Layer  Model  Input  for  2  November  1999,  French  Alps 


base  (km) 

-A (s'1) 

ddd/ff  f  /.  m/s) 

Layer  1 

2.5 

0.012 

220/15 

Layer  2 

7.0 

0.009 

220/21 

Layer  3 

12.0 

0.020 

220/31 

Note  that  the  wind  direction  is  constant  with  height  while  the  speed  varies  between 
layers. 

Figure  6  depicts  the  model  height  perturbation  solution  r/(x,y )  in  meters  at  z  =  7.5  km 
MSL  in  the  horizontal  plan  view  (top)  and  q  (x,y,z)  along  a  SW  -  to  -  NE  cross-section 

(bottom).  The  cross-section  horizontal  position  is  depicted  by  the  blue  dashed  line  on  the 
plan  view.  The  axis  of  the  bow  wave  pattern  resulting  from  the  terrain  height  distribution 
lies  parallel  to  the  mean  wind  flow  and  a  variety  of  transverse  wave  patterns  are  also 
evident  in  the  flow.  The  vertical  cross-section  in  the  lower  half  of  the  figure  indicates  the 
amplitude  of  the  height  perturbations  at  each  level  and  the  tendency  for  these 
perturbations  to  tilt  upstream  (to  the  left  in  the  bottom  figure). 


5.3.3  COMPARISON  WITH  BALLOON  MEASUREMENTS 

To  gain  further  appreciation  into  the  validity  of  the  Fourier  calculations,  a  comparison  of 
model  vertical  velocity  w  with  balloon-derived  air  vertical  velocity  was  accomplished. 
Basic  states  for  the  Fourier  calculations  were  constructed  using  Global  Forecast  Spectral 
Model  archives  for  two  different  times  and  locations  where  high  resolution  ( —  /1 .5  second 
sample  rate)  balloon-borne  rawinsonde  measurements  were  made  in  the  vicinity  of 
mountains.  The  locations  and  times  for  the  balloon  launches  were  (1)  the  island  of 
Hawaii  with  balloon  launch  from  Bradshaw  Army  Airfield  (19°47’N,  155°33'  W)  on 
12  December  2002  near  06  UTC  and  (2)  the  Provence  Region  of  France  in  the  southern 
Alps  with  balloon  launch  from  the  OHP  Observatory  (44°  N,  5°  42'  E)  on  24  November 
2004  around  00  UTC. 


km  (north) 


Az  at  7.5  km  MSL 
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Figure  6.  Height  perturbations  //(x,y)  displayed  for  French  Alps  MAP  case  with  basic  state 

as  described  in  Table  1  for  z  =  7.5  km  MSL  (top)  and  for  SW-  NE  vertical  cross-section  at 
bottom.  The  cross  section  is  along  the  line  denoted  by  the  dashed  blue  line  on  the  plan  view. 


The  basic  state  atmosphere  over  Hawaii  was  characterized  by  weak  northeast  trade  flow 
at  low  altitudes  backing  toward  stronger  northwest  flow  above  10  km.  A  trade  inversion 
existed  with  the  inversion  base  near  1.5  km  MSL.  For  the  Provence  atmospheric 
background,  steady  northerly  flow  prevailed  at  all  levels  with  some  weak  backing  of  the 
wind  with  height.  The  basic  states  for  the  two  locations  are  listed  in  Tables  3  and  4 
where  U  represents  the  west-to-east  wind  speed  component  and  V  represents  the  south- 
to-north  wind  speed  basic  state  component. 

Table  3.  Three-Layer  Model  Specifications  for  Hawaii 
12  December  2002  06  UTC 


site  elevation:  1887  m 

>■55 

II 

O 

o 

(boundary  reflectivity) 

;  stable  layer : 

0  -  1  km  MSL 

base  (km) 

Kish 

U  (m/s) 

V  (m/s) 

Layer  1 

1.0 

0.010 

-  10.0 

-5.0 

Layer  2 

9.0 

0.008 

+20.0 

-10.0 

Layer  3 

16.0 

0.020 

+  2.0 

-2.0 

Table  4. 

Three-Layer  Model  Specifications  for  Provence 

24  November  2004  00  UTC 

site  elevation:  650  m 

o 

as 

o' 

II 

o 

(boundary  reflectivity) 

;  stable  layer : 

0  -  0.6  km  MSL 

base  (km) 

U  (m/s) 

V  (m/s) 

Layer  1 

0.6 

0.016 

-0.7 

-3.9 

Layer  2 

3.1 

0.008 

+  1.4 

-7.8 

Layer  3 

10.5 

0.020 

+  1.6 

-9.1 

A  high-pass  filter  technique  was  used  to  calculate  the  observed  vertical  velocity 
perturbations  about  the  mean  balloon  rise  rate  (Murphy,  2006)  using  standard  digital 
processing  techniques.  The  Nyquist  frequency  of  the  balloon  sampling  rate  is  much 
higher  than  the  spatial  frequencies  associated  with  the  mountain  waves.  The  mean 
balloon  rise  rate  is  determined,  using  the  hydrostatic  assumption,  from  the  sonde  pressure 
time  tendency  Ap/At .  Air  pressure  and  temperature  are  sampled  at  1-  to  2-second 
intervals  during  balloon  ascent.  Horizontal  wind  velocities,  reported  at  a  comparable 
sampling  interval  by  the  rawinsonde  using  GPS-derived  Doppler  shift  and  smoothed 
using  a  Gaussian  window  filter,  are  also  used  to  determine  the  three-dimensional  balloon 
trajectories.  The  comparison  between  balloon  and  model  vertical  velocities  was 
performed  from  the  surface  to  15  km  above  sea  level  (MSL).  For  these  cases,  the 
horizontal  balloon  drift  during  the  ascent  to  15  km  was  on  the  order  of  20  to  40  km. 

The  model  vertical  velocity  w  along  the  balloon  trajectory  s(x,y,z)  used  in  the 

comparison  is  interpolated  from  the  three-dimensional  Fourier  model  output  grid  using  a 
vertical  grid  spacing  Az=  50  m  and  horizontal  grid  spacing  Ax  =  Ay  =  1  km.  15  km  MSL 
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Figure  7a.  Terrain  relief  overlain  with  balloon  trajectories  tracked  to  15  km  for  Hawaii,  12 
December  2002. 
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Figure  7b.  3  Layer  Model  vertical  velocity  contours  over  Hawaii  at  12  km  MSL.  The 
contour  interval  is  0.4  m/s.  The  red  ring  marks  the  balloon  launch  location. 
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Figure  8a.  Terrain  relief  overlain  with  balloon  trajectories  tracked  to  15  km  in  France,  24  November 
2004. 


Figure  8b.  3  Layer  Model  vertical  velocity  contours  over  Provence  at  12  km  MSL.  The 
contour  interval  is  0.5  m/s.  The  red  ring  marks  the  balloon  launch  location. 


hhh02001  Vertical  Velocity  compared  to  Model 


Figure  9a.  Balloon-derived  vertical  velocity  (blue)  compared  with  3  Layer  model  output 
(red)  for  sonde  launch  at  12  December  2002  0500  UTC  from  Bradshaw  AAF,  HI 


hhh02002  Vertical  Velocity  compared  to  Model 


Figure  9b.  Balloon-derived  vertical  velocity  (blue)  compared  with  3  Layer  model  output 
(red)  for  sonde  launch  at  12  December  2002  0657  UTC  from  Bradshaw  AAF,  HI 
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Figure  10a.  Balloon-derived  vertical  velocity  (blue)  compared  with  3  Layer  model  output 
(red)  for  sonde  launch  at  23  November  2004  2335  UTC  from  OHP  Observatory,  France 
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Figure  10b.  Balloon-derived  vertical  velocity  (blue)  compared  with  3  Layer  model  output 
(red)  for  sonde  launch  at  24  November  2004  0406  UTC  from  OF1P  Observatory,  France 


is  set  as  the  top  level  of  model  calculation  for  these  comparisons,  which  equates  to  a  less 
than  50  minute  period  using  a  mean  balloon  ascent  rate  of  5  m/s. 


Figures  7  and  8  depict  the  terrain  relief,  horizontal  balloon  trajectories,  and  model 
vertical  velocity  output  at  12  km  MSL  for  Hawaii  and  Provence,  respectively.  The  same 
axis  scale  reference  is  used  for  the  terrain  and  vertical  velocity  plots  in  each  figure. 

These  figures  illustrate  the  complexity  of  both  the  underlying  terrain  and  the  associated 
model  vertical  velocity  fields  generated  using  the  terrain  and  atmospheric  background. 

Figures  9  and  10  show  the  comparison  of  the  model  and  balloon-derived  vertical 
velocities  for  the  Hawaii  and  Provence  launches,  respectively.  For  the  Hawaii  case 
(Figure  9),  the  model  magnitudes  were  much  larger  than  the  balloon-derived  values 
below  8  km  but  showed  better  agreement  in  sign  and  magnitude  between  9  and  15  km 
MSL.  Below  15  km,  the  maximum  measured  vertical  velocities  are  on  the  order  of  +/-  1 
m/s  while  the  model  predicted  velocity  magnitudes  are  on  the  order  of  +/-  2  m/s.  There 
was  also  better  agreement  in  the  sign  of  the  vertical  velocity  above  8  km  compared  to 
lower  altitudes. 

For  the  Provence  case  study,  reasonable  agreement  in  magnitude  and  sign  were  found 
between  the  observed  and  model  vertical  velocities  (Figure  10)  although  the  model  was 
somewhat  out-of-phase  with  the  observed  during  the  first  launch  (2335  UTC).  For  the 
two  launches,  the  magnitude  of  the  observed  and  model  vertical  velocities  varied  between 
approximately  +/-  1  m/s.  Due  to  poor  wind  data  acquisition  during  the  second  launch 
(0406  UTC),  the  balloon  trajectory,  and  thus,  the  model  comparison,  ends  just  above  5.4 
km  MSL. 


6.  BROUTMAN  CONTINUOUS  MOUNTAIN  WAVE  FORECAST 
MODEL  3  (MWFM-3) 

Three  principal  journal  articles  describe  MWFM-3,  the  intended  successor  to  MWFM-2 
(Marks  and  Eckermann,  1995).  They  are  Broutman  et  al,  2002  ;  Broutman  et  al,  2003  ; 
and  Broutman  et  al,  2006.  The  first  article  (2002)  deals  with  stationary,  hydrostatic, 
vertically-propagating  wave  modes  while  the  second  article  (2003)  deals  primarily  with 
trapped  non-hydrostatic  modes  involving  steady-state  solutions.  The  most  recent  article 
(2006)  describes  transient  or  time-dependent  solutions.  These  solutions  are  all 
demonstrated  using  a  shallow  quasi-Gaussian  shaped-hill  with  a  height  of  100  m  for  the 
lower  boundary.  This  shape,  sometimes  referred  to  as  the  “Witch-of- Agnesi”  and  easily 
represented  in  the  wave  number  domain  by  an  algebraic  expression,  is  used  to  induce 
small  magnitude  stationary  perturbations  in  the  height  field.  From  these  height 
perturbations,  vertical  and  horizontal  velocities  are  also  diagnosed  using  analytical 
expressions. 


6.1  Governing  Equations  for  MWFM-3  (after  Broutman  et  al,  2003): 


The  intrinsic  frequency  a  and  dispersion  relation  for  the  vertical  wave  number  777, 
respectively,  are: 


a(kx,ky,z)  =  -kxU(z)-  ky  V  (z) 


(k2x+k^N{z) 
kl  +k2y  +m(kx,ky,zj 


(rad/s) 


For  non-hydrostatic  wave  modes: 


(13) 


777 


(kx,ky,z)=  [k2x  +  k2y ) 


,N2  -<j2 


(rad/m)  (14) 


For  hydrostatic  wave  modes: 


m(kx,ky,z)=  {k;+k2) 


N2 


(rad/m)  (15) 


The  height  perturbation  //  is  composed  of  the  sum  of  vertically  propagating  fjpr  and 
trapped  wave  fjtr  components  over  all  wave  numbers: 

ij{kx,ky,z)  =  i]pr  (, kx,ky,z )  +  ijtr  (kx,ky,z)  (m3/rad) 


From  this,  the  vertical  and  horizontal  velocity  in  the  transfonn  domain  can  then  be 
determined 


w(kx,ky,z)  =  -icTTj(kx,ky,z) 
-ikjnw 


(m3  s'1) 


i(kx,ky,z)  = 


(k>K) 


The  buoyancy  frequency  turning  point  height  zt  =  f  {k,  N,  U)  is  a  function  of  horizontal 
wave  number  or  mode,  static  stability,  and  horizontal  wind  velocity  and  is  defined  locally 
to  be  where  intrinsic  frequency  is  equivalent  to  the  static  stability  or  a-N.  This  also 
serves  as  a  caustic  point  (Lighthill,  1978)  since  the  turning  point  acts  as  the  boundary  for 
a  trapped-wave  region. 


According  to  Lighthill  (1978),  a  caustic  is  a  boundary  between  a  region  with  a 
complicated  wave  pattern  due  to  interference  between  two  groups  of  groups  of  waves, 
and  a  neighboring  region  including  no  waves.  Both  critical  layers  and  turning  points  act 
as  caustics,  although  turning  points  only  occur  with  non-hydrostatic  modes.  The 
MWFM-3  model  described  here  can  only  handle  the  cases  of  zero  or  one  turning  point  in 
the  vertical  for  each  wave  number  mode. 


6.2  Vertically  propagating  wave  mode: 


The  vertically-propagating  wave  mode  is  defined  as  the  contribution  to  the  upward 
propagating  wave  field  for  each  mode  at  or  below  a  turning  point  height  zt.  (z,  does  not 
necessarily  exist  for  each  mode,  e.g.  for  hydrostatic  modes.)  The  phase  angle  for 
vertically  propagating  waves  is 

<t>p  ( K , ky )  =  - Jo‘  m ( k,  ,ky,z)dz  (radians)  ( 1 6) 

where  the  integral’s  sign  is  assigned  the  following  meanings: 


(-)  sign  for  m  indicates  upward  group  velocity  and  downward  phase  velocity 
(+)  sign  for  m  indicates  downward  group  velocity  and  upward  phase  velocity 


Differentiating  Eq.  (13)  with  respect  to  m,  the  vertical  wave  group  velocity  cg  is  written 
as: 


cg(kx,ky,z)  = 


da  ~{kl  +  k;  f  ■  N (z)  •  m 


(m  s'1) 


dm  (kl+K+mT 

The  vertical  flux  of  wave  action  G(kx,ky,z )  is  expressed  in  terms  of  static  stability 

square  of  the  Brunt-Vaisala  frequency,  N2 ),  vertical  group  velocity,  and  intrinsic 
frequency: 


(17) 

(the 


G(kx,ky,z)  =  N2(z)cg3/a(kx,ky)  (kg  in'2  s'2) 


(18) 


The  expression  for  height  perturbation  for  vertically  propagating  modes,  using  the 
general  form  of  Eq.  (1),  then  becomes 


dpr(kx,ky,z)  =  h(kx,ky)[GjG'fl  e4p  (m3  rad'1)  (19) 

where  h{kx,ky)j  =  ^—^  J  h{x,y)dk*xek>y  dxdy  (m3  rad'1)  and  the  “0”  subscript 

denotes  the  expression  for  a  variable  that  is  valid  at  z  =  0. 


For  vertically-propagating  wave  modes,  solution  methods  are  straightforward  and  involve 
simply  the  vertical  integration  of  the  dispersion  relation  m  from  the  surface  to  the  desired 
altitude  z  .  Examples  of  these  solutions  are  displayed  in  Figure  11.  A  20  km  radius  hill, 
centered  at  the  origin  with  a  height  of  100  m,  is  used  for  the  lower  boundary.  Contour 
plots  are  produced  for  the  two-dimensional  distribution  of  vertical  and  horizontal  velocity 
at  altitudes  equal  to  2. 1  km  and  6.3  km,  and  compared  with  published  results  from 
Broutman  et  al  (2002).  The  blue  values  in  parentheses  represent  the  published  maxima 
and  minima  while  the  black  values  represent  the  AFRL  model  calculation  results.  The 
phase  functions  used  in  these  examples  are  evaluated  analytically.  Close  agreement 
between  published  and  calculated  magnitudes  show  the  similarity  in  application  of 
Fourier  methods  even  though  the  grid  dimensions  differed  significantly  —  28  for 
Broutman  et  al  (2002)  vs.  2 10  for  the  results  displayed  in  Figure  11. 
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Figure  11.  Calculations  of  hydrostatic  vertically-propagating  modes  using  MWFM-3 
Values  in  parentheses  are  published  results  (Broutman  et  al,  2002) 

numerical  specifications: 

Ax  =  Ay  =  1000  m 
bandwidth:  6.28  x  10"3  rad/m 
wave  number  grid:  1024  x  1024 


atmospheric  specification: 

h(z)  -U0  cos (nz/H) ;  v(z)  =E/0  sin (nz/H)  ;  H  =  12  km  ;  U0=  10  m/s 
N  =  0.01 13  s"1  ;  Richardson  number  ~  19  ;  hill  radius  =  20  km,  height  =  100  m 


(d) 


6.3  Vertically-trapped  waves:  fjtr 


For  vertically-trapped  waves,  the  “long-time”  steady-state  solution  requires  invocation  of 
the  Airy  function  Ai(p)  where  the  phase  functions  form  its  argument,  p  ,  and  which  may 
be  complex.  The  Airy  function  and  its  determination  are  described  in  Section  6.4  and  in 
Appendix  B.  Solutions  which  correct  for  the  trapped  waves  or  caustics  are  referred  to  as 
“uniform”  solutions  (Broutman  and  Rottman,  2004)  due  to  their  unifonnly  valid  behavior 
in  the  vicinity  of  singularities. 

In  a  general  form,  the  trapped  wave  solution  is: 

nlr(kx,ky,z)  =  h(kx,ky,z)[GjGf  [p/pQf  x  (20) 

The  phase  functions  <j>  are  expressed  as 

(j\  =  -J  '  m  dz  for  a  ray  incident  upon  the  caustic  level  z, 

(f>2  =  ±J  m  dz  for  a  ray  reflected  from  the  caustic  level  zt 

where 

p  =  -$>)]^  >  Ai(p)  is  the  Airy  function,  and  p  can  be  real  or  complex 

and  p  =  0  if  z  =  zt 
p  <  0  if  z  <zt 
p>  0  if  z>  z, 

r  3 

p  =  p0 for z  =  0,  p0  =  -  -<t>x 

For  the  solutions  described  below,  p  is  treated  as  a  real  function  and  p0  is  complex. 

The  phase  integral  sign  convention  follows  the  same  relation  as  was  used  for  vertically- 
propagating  modes: 

+  upward  phase  velocity 
downward  group  velocity 

-  downward  phase  velocity 
upward  group  velocity 

The  relevant  features  of  Broutman  et  al  (2003)  for  stationary  trapped  wave  modes  were 
replicated  to  gain  understanding  into  the  trapped  wave  mode  behavior.  The  calculations 
use  the  “Case  II”  specification  of  Wurtele  et  al  (1987).  This  case  uses  a  linear  increase  of 
zonal  wind  speed  with  height  and  a  fixed  value  for  static  stability  (N  =  0.01  s'1).  A  2.5 
km-radius  hill  centered  at  the  origin  with  height  of  100  m  is  used  for  the  lower  boundary. 
All  calculations  are  for  trapped  wave  modes  and  use  complex  functions  except  where 
indicated.  The  essential  results  are  presented  in  Figures  12  through  16. 


Figure  12  shows  that  calculations  of  turning  point  height  zt  as  a  function  of  horizontal 
wavelength  compare  well  with  the  published  results.  For  calculations  at  all  wavelengths, 
if  zt  =  0,  then  the  eigenfunction  77,,.  is  treated  as  zero  at  that  wavenumber.  The  turning 
point  height  is  also  treated  as  zero  in  MWFM-3  once  it  becomes  less  than  zero.  The 
turning  point  height  will  eventually  become  zero  or  negative  when  the  spatial  wave 
number  exceeds  a  certain  value,  as  shown  in  Figure  12.  For  the  MWFM-3  model 
calculations,  these  zero  or  negative  turning  point  height  wavenumber  values  are  skipped. 

The  nonnalized  vertical  eigenfunction  calculations  for  trapped  wave  modes  using 
MWFM-3  plotted  in  Figures  13a  and  13c  agree  reasonably  well  with  those  published. 
However,  there  is  a  sign  reversal  in  Figure  13b  which  was  not  replicated  and  which  may 
contribute  to  some  of  the  differences  seen  in  Figure  14.  The  eigenfunctions’  appearance 
also  became  slightly  smoother  with  altitude  when  the  phase  function  integrals  were 
evaluated  analytically  rather  than  numerically. 

In  Figure  14,  the  eigenfunction  modulus  contour  patterns,  as  a  function  of  spatial 
wavenumber,  match  but  there  is  slightly  more  damping  of  the  peaks  in  the  calculations 
for  the  uniform  solution  compared  to  published  values.  Here,  the  15  km  mode  magnitude 
dominates  the  33  km  mode.  The  uniform  solution  damping  with  height  is  due  to  the 
addition  of  a  small  imaginary  wavenumber  component,  ki  =  5  x  10'6  rad/m  .  Additional 
experiments  (not  shown)  indicate  the  eigenfunction  magnitude  and  wavenumber 
positions  are  somewhat  sensitive  to  the  size  of  the  imaginary  wavenumber  component. 

Figure  15  shows  that  the  vertical  velocity  cross-sections  for  the  1-D  Fourier  transfonn 
solution  match  well  the  published  values.  The  maximum  published  values  are  shown 
parenthetically.  The  vertical  velocity  is  calculated  at  1  km  vertical  increments.  The 
calculated  vertical  cross-section  using  the  1-D  Fourier  transform  and  the  2-D  inverse 
Fourier  transform  displayed  in  plan  fonn  (Figure  16  a)  matches  the  dominant  ~  15  km 
mode  shown  in  Figure  14  and  Wurtele  et  al’s  linear  analysis  for  case  II.  Evidence  of  a  30 
-35  km  mode  is  seen  near  an  altitude  of  15  km  in  Figure  15. 

Figure  16  shows  a  comparison  of  the  MWFM-3  trapped  wave  vertical  velocity  pattern  for 
the  2-D  plan  view.  Broutman  et  al  (2003)  attribute  the  ship  bow  behavior  and  modulus 
diminishment  with  distance  downstream  of  the  hill  to  the  imaginary  wave  number 
component.  Both  calculations  show  excellent  agreement  in  amplitude  and  phase.  The 
published  maximum  and  minimum  values  are  shown  in  parentheses. 

Two  modifications  were  required  beyond  what  is  cited  in  the  Broutman  et  al  (2003) 
article  to  produce  the  trapped  wave  modes  that  dissipate  downstream  with  distance. 

These  were  detennined  from  correspondence  with  Dr.  Broutman: 

The  primary  modification  was  that  the  phase  function  po  should  be  maintained  as  a 
complex  argument  of  the  complex  Airy  function  Ai(p0)  in  Eq.  (20).  To  avoid 
singularities  associated  with  the  resonant  modes  or  zeroes  of  the  Airy  function,  a  small 
imaginary  component  0.01/  is  added  to  po.  However,  p  is  maintained  as  the  real 
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Figure  12.  Turning  point  height  zt  as  a  function  of  wave  number. 

Atmospheric  specification:  (Case  II  from  Wurtele  et  al,  1987) 
U=U0-(l+z/L}  TV  =  0.01s-1  a  =  2500m  li  =  1 00 m 

U0  =  10  ms  1  ;  L  —  4  km  Richardson  number  =  16 
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Figure  13a.  Comparison  of  vertical  eigenfunction  Rcy]u/  hj  for  A.  =10  km  as  a  function  of  height 
using  MWFM-3  (2-D)  for  trapped  wave  modes.  Dashed  line  marks  turning  point  height  zt. 

Lower  boundary  specification:  h(k)  =  h0ae  a'k'  ht]  —  1 00  m  (Witch-of-Agnesi) 

Published  solution  (Broutman  et  al,  2003)  is  at  top. 
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Figure  13b.  As  in  figure  13a,  except  Rc(  i]N  / h  j  for  1=20  km.  Vertical 
dashed  line  indicates  turning  point  height  zt. 
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Figure  13c.  As  in  figure  13a  but  for  comparison  of  vertical  eigenfunction 
Re  ( fju  fh  I  for  /.=  40  km.  Vertical  dashed  line  indicates  turning  point  height  zt. 
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Figure  14.  From  Broutman  et  al  (2003)  representing  uniform  2-D  MWFM-3  solution 

(m2  rad"1)  for  p  =  —  1.  Parentheses  denote  eigenfunction  peak  magnitude  (m2  / 

rad12).  Atmospheric  and  terrain  specifications  used  are  as  listed  in  Figures  12  and  13. 
Published  solution  is  at  top. 
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Figure  15.  Comparison  of  two-dimension  cross-section  for 
MWFM-3  steady-state  version  with  trapped  waves. 

ki=  5  x  10'6  rad/m  wmax=  +  0.39  m/s  (+0.38) 

s  =  1024  km  wmin  =  -  0.39  m/s  (-0.37) 

n=  1024  Aw  =  0.04  m/s 

Az  =  1  km 
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Figure  16.  MWFM-3  trapped  wave  solutions  Published  value  at  top 
ki  =  5  x  10'6  for  kx  only  same  terrain  and  atmospheric  specification  as  Fig.  13 
s  =  1024  km  z  =  2.5  km 

n  =  1024  x  1024  wmax  =  +  0.25  m/s  (+0.21) 

wmin=  -  0.24  m/s  (-0.18) 

Aw  =  0.02  m/s 


argument  of  the  Airy  function.  A  secondary  modification  involved  evaluation  of  the 
phase  integrals  analytically  instead  of  numerically.  The  dispersion  relation  for  this  case 
is  in  a  form  which  permits  ready  evaluation  using  integral  tables.  Once  these 
modifications  were  accomplished,  the  results  in  Figures  15  and  16  matched  the  published 
results.  This  behavior  suggests  the  MWFM-3  results  for  trapped  wave  modes  can  be 
sensitive  to  phase  function  evaluation. 

The  trapped  mode  results  displayed  in  Figures  15  and  16  can  also  be  tested  and 
reproduced  using  a  non-linear  numerical  model  given  the  same  initial  conditions.  In 
Broutman  et  al  (2003),  a  nonlinear,  Boussinesq,  second-order  finite  difference  numerical 
model  is  run  out  to  1.5  hours  until  a  steady-state  solution  is  obtained.  It  uses  a  “sponge 
layer”  for  x  >  80  km  and  z  >  35  km.  The  results  from  the  nonlinear  model  and  linear 
model  match  almost  exactly. 

6.4  Use  of  the  Airy  Function  in  MWFM-3 

MWFM-3  solutions  for  trapped  wave  modes  require  evaluation  of  the  Airy  function,  also 
referred  to  as  the  Airy  integral,  for  real  and/or  complex  arguments.  The  Airy  integral  can 
be  applied  in  the  vicinity  of  caustics  to  treat  the  breakdown  of  ray  wave  theory  (Lighthill, 
1978).  Its  analytical  evaluation  procedures  are  described  in  Appendix  B.  The  Airy 
function  Ai(z)  can  be  calculated  very  accurately  using  ascending  series  where  z  =  x  +  iy  is 
complex.  However,  the  series’  relatively  slow  convergence  dictates  a  much  quicker 
calculation  method  when  potentially  millions  of  evaluations  must  be  repeatedly  made. 

For  parameterization  of  the  Airy  function  involving  only  the  real  argument  x,  a  three- 
zoned  approximation  can  be  used  where  (1)  auxiliary  functions  are  employed  toward  the 
extrema,  (2)  an  approximation  from  Lighthill  (1978)  is  used  for  significant  negative  x, 
and  (3)  a  polynomial  fit  is  used  in  the  vicinity  of  x  =  0.  Both  the  ascending  series 
solution  and  the  fast  parameterization  are  plotted  in  the  first  figure  of  Appendix  B  to 
illustrate  their  equivalence. 

A  machine-oriented  algorithm  (Amos,  1986)  is  used  for  fast  calculation  of  the  Airy 
function  with  complex  arguments.  Its  evaluation  in  the  complex  plane  compares  well 
with  those  obtained  using  the  series  expansion,  as  displayed  in  the  modulus  and  phase 
contour  charts  of  figures  B2  and  B3  in  Appendix  B. 

6.5  Numerical  Evaluation  of  the  Phase  Integral  for  MWFM-3: 

Implementation  of  MWFM-3  for  use  with  routine  weather  data  requires  the  capability  to 
numerically  evaluate  the  phase  integrals.  Figures  17  and  18  depict  plots  of  vertical 
velocity  (m/s)  calculated  at  6.3  km  using  the  hydrostatic  assumption  for  the  dispersion 
relation,  constant  Brunt-Vaisala  frequency,  and  the  horizontal  wind  velocity  for  two  cases 
—  one  where  the  wind  vector  is  backing  with  height  and  one  where  the  zonal  wind  is 
constant  with  height.  The  lower  boundary  terrain  is  a  200  meter  bell-shaped  hill  with  a 
half- width  of  20  km  (see  figure  19).  The  phase  integral  is  calculated  numerically  over  64 


levels  (A z  =100  m)  using  the  trapezoidal  rule  on  a  1024  x  1024  horizontal  mesh  (A  ~  1 
km)  in  lieu  of  an  analytical  evaluation  of  the  phase  integral.  Critical  levels,  where  the 
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z  =  6.3  km  MSL,  Az  =  100  m,  64  levels 
max  vertical  velocity  =  0.108  m/s 
min  vertical  velocity  =  -  0.066  m/s 
hydrostatic  dispersion,  wind  backing  with  height 
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Shutts  and  Gadian,  1999 
(analytical)  (+  0.08,  -0.07) 


Doyle  and  Jiang,  2006 
(analytical)  (+  0.09,  -0.06) 


Figure  17.  Plot  of  hydrostatic  mountain  wave  solutions  with  base  state  wind  backing  with  height  as 
specified  above.  Dashed  line  on  the  Shutts  and  Gadian  plot  shows  the  wind  direction  at  z  =  6.3  km. 
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Figure  18.  Plot  of  hydrostatic  mountain  wave  solutions  with  constant  zonal  wind  with  height  as 
specified. 


intrinsic  frequency  is  zero,  are  skipped.  1004  critical  levels  were  detected  for  the  backing 
wind  case  out  of  a  total  ~  67  x  106  points. 

The  examples  depicted  in  Figures  17  and  18  were  originally  published  by  Shutts  and 
Gadian  (1999)  but  other  authors  have  replicated  these  examples  including  Tanushev  et  al 
(2007)  using  the  Gaussian  beam  method  and  Doyle  and  Quing  (2006)  using  a  modified 
version  of  Smith’s  model.  At  the  bottom  of  the  figures,  the  calculated  results  for  wind 
backing  with  height  are  compared  with  their  results  as  well  as  Shutts  and  Gadian’ s  non¬ 
hydrostatic  numerical  simulations  and  hydrostatic  analytical  solutions.  The  maximum 
and  minimum  vertical  velocities  obtained  by  the  different  investigators  generally  agree. 


x  (km) 

Fig.  19  Bell-shaped  terrain  used  in  model  calculations  (meters),  for  figures  17  and  18,  with 
200  meter  height  and  20  km  half-width.  Inner  contour  intervals  are  approximately  20  meters. 


6.6  MWFM-3:  Discrete,  Steady-State,  Linear,  Hydrostatic  and  Non¬ 
hydrostatic 

Broutman’s  linear  steady-state  MWFM-3  model  for  hydrostatic  and  non-hydrostatic 
solutions  can  be  expressed  in  fully  discrete  form.  The  distinguishing  features  between 
this  and  results  presented  above  are  the  provisions  for  (1)  discrete,  multi-level  input  for 
the  basic  state,  such  as  provided  by  rawinsonde  or  numerical  models,  including  the 
digitized  terrain  data  base,  and  (2)  a  switch  for  calculating  hydrostatic  or  non-hydrostatic 
solutions  including  trapped  wave  modes. 


wmax  = +.102383,  wmi„  = -0.062303  z  =  6.3  km 
hm  =  200  m  ;  Gaussian  hill  !4  width  =  20  km 
N  =  0.01  s'1,  U=  10  m/s 
Vsfc  =  0,  dv/dz  =  3  x  10'3  s'1 


Shutts  and  Gadian,  1 999 
(analytical)  (+  0.08,  -0.07) 


Figure  20.  Hydrostatic  (wind  backing  with  height)  .  Fully  Discrete  version 


wmax  =  +0.25  m/s  ,  wmin  =  -  0.24  m/s 

contour  interval  =  0.02  m/s  hmax  =  100  m  Arcs  ~  15  km 

Gaussian  half-width  =  2.5  km  ,  z  =  2.5  km 

sfc:  Uo  =  10  m/s  dU/dz  =  10/4000  s'1  N  =  0.01  s'1 


wmaz=  +0.21  m/s,  wm;n=  -0.18  m/s  ,  A^es  ~  15  km 

analytical  inputs  for  shear,  stability,  turning  points  (Broutman  et  al,  2003) 


wmax  =  +0.25  m/s,  wmin=  -0.21  m/s  (Broutman  et  al,  2003) 
non-linear,  Boussinesq,  free-slip  lower  boundary,  sponge  layer  for  x  >  80  km 


Figure  22.  Non-hydrostatic  (constant  wind  shear  with  height;  fixed  static  stability)  Fully  discrete 


Analytical  functions  and  dependencies  involving  the  basic  state  (wind  velocity  and 
temperature)  and  phase  behavior  are  removed.  A  two-dimensional  lower  boundary  is 
defined  in  real  space  and  then  transfonned.  Phase  integrals  are  evaluated  numerically 
using  the  trapezoidal  rule.  The  existence  of  turning  points  and  critical  levels  are  checked 
for  and  determined  objectively  at  each  level  for  all  wave  numbers.  An  upper  altitude 
limit  of  30  km  is  imposed  for  turning  point  checks,  consistent  with  the  upper  limit  of 
typical  NWP  and  rawinsonde  data.  Altitude  increments  can  vary  within  the  vertical 
column.  For  the  examples  here,  a  uniform  height  increment  of  100  m  is  used. 

These  particular  examples  employ  numerically-determined  turning  points  and 
numerically-evaluated  phase  integrals  for  hydrostatic  (Figures  20  and  21)  and  non¬ 
hydrostatic  trapped  wave  test  case  (Figure  22)  using  Broutman’s  steady-state  model 
(Broutman  et  al,  2003).  For  evaluation  purposes,  a  turning  point  exists  when  the  real 
part  of  the  vertical  wave  number  (dispersion  relation)  changes  sign  between  levels,  or  its 
modulus  is  zero  at  a  level.  Critical  levels,  where  the  intrinsic  frequency  is  zero,  and  “zero 
altitude”  turning  points,  are  skipped  for  the  Fourier  transfonn  evaluation. 

Comparisons  of  MWFM-3  steady-state  discrete  solutions  (in  the  upper  figure  half)  with 
published  results  (in  the  lower  half)  are  shown  for  the  hydrostatic  and  non-hydrostatic 
cases.  Solutions  for  vertical  velocity  are  calculated  on  a  1024  x  1024  horizontal  mesh 
using  a  mesh  spacing  of  1  km.  Wind  profiles  are  treated  effectively  as  continuous 
functions  using  constant  values  of  speed  or  shear  with  height.  Stability  values  can  also  be 
allowed  to  change  abruptly  between  levels  for  the  non-hydrostatic  model  comparisons. 
However,  such  a  stability  profile  behavior  may  not  always  be  appropriate  for  a  multi¬ 
level  discrete  approach. 

One  deficiency  in  the  MWFM-3  formulation  as  published  is  the  neglect  of  boundary 
layer  wave  absorption  associated  with  the  surface  stagnant  layer,  as  in  the  manner 
discussed  by  Smith  et  al  (2002). 


7.  TRANSIENT  WAVE  SOLUTIONS  (MWFM-3): 

(Broutman  et  al,  2006) 

The  general  transient  or  time-dependent  trapped  wave  solution  for  fjtr ,  described  in 
Broutman  et  al  (2006),  is: 

W=7„  +{M, -M2)h[GjGfe^+M^  =  7jtr(kx,ky,Z,  t)  (m3rad'2)  (21) 

where  fju  =  2z';r^(-r)^ A- Ai(r)e^~^SMi  (22) 

A  =  h  [G0/G]':  (wave  action  density),  Mi  and  M2  are  the  trapped  wave  counters 
G  =  pN2cg  la  (kg  m"2  s'2  rad'1),  p  here  being  air  density  (kg/m3) 

Ai(r)  is  the  Airy  function  for  the  argument  r 

h{kx,ky )  represents  the  topography  of  the  lower  boundary  in  the  wave  number 
domain.  For  a  bell-shaped  hill,  h(kx,ky)  =  h0a2e~a{ICx+ky)/2 


1  _  &iM 2* 

— — —  =  0  if  M 2  =  0  is  the  summation  expression  for  trapped  waves.  The 


time-dependent  trapped  wave  behavior  and  related  nomenclature  are  discussed  further 
below.  The  group  velocity  is,  as  described  earlier  in  Eq.  (17), 

Cg=^m=  +  k>' ^  ' m(z') '/(*’  +  ky  +  m (z)2 T  ’ 

For  the  lower  boundary  condition,  the  height  field  perturbation  is  set  equivalent  to  the 
topography,  or  ij {kx,ky,0)  =  h[kx,ky\.  The  Airy  function  argument  r  is  complex.  In  the 


transient  model  implementation,  however,  both  Ai(r)  and  its  argument  r  are  treated  as 
real  functions. 


The  phase  relations  for  the  transient  wave  solution  are  similar  to  those  for  the  steady-state 
solution.  For  the  vertically  propagating  waves, 

</>(kx,ky,z}  =  -^  m(kx,ky,z}dz  (radians)  (23) 

For  the  trapped  waves, 

let  Q>  =  2<j>-n£,  where  (j>  =  -^mdz  ;  and  ^  =  \[cj)l+ r  =  -[f(^l -^2)]^ 

If  z,  <  z  where  z  is  the  level  of  analysis,  then  the  expression  for  r  will  have  a  positive 
sign.  If  zt-z,  then  r  -  0. 


Let  tj)2  represent  the  phase  for  waves  reflected  from  the  caustic  and  ^  represent  the  phase 
for  waves  incident  upon  the  caustic.  These  phase  functions  are  defined  as: 

(f>2  -  cj)  -  J  '  m  dz  -  -  J ^  m  dz  -  J  '  m  dz  (radians),  the  phase  for  waves  reflected  from  the 
caustic,  and 


(j>x  =  -J"  m  dz  ,  the  phase  for  waves  incident  upon  or  below  the  caustic 

z  is  the  analysis  level,  and  z,  is  the  turning  point  height  (m),  that  is,  the  height  where  the 
vertical  wave  number  m  becomes  less  than  or  equal  to  zero. 


The  relevant  wave  dispersion  relations  for  the  transient  trapped  solutions  are  the  same  as 
the  non-hydro  static  steady-state  formulation: 


m  -  k. 


and 


N2-ct2 


(rad/m)  ;  a  =  -kJJ  -  kyV  (rads’1)  ;  kh  =  [k2  +  k2v  (radnT1)  (24) 


h{kx,ky^-^—^  |  h ( x,y)e‘k'xe,k’ydx dy  (m3  rad’1) 


(25) 


For  the  example  described  below  in  section  7.2,  U (z)  =  U0  (l  +  z/L )  (m  s’1) 


According  to  Broutman  et  al  (2006),  no  imaginary  component  need  be  added  to  the 
horizontal  wave  number  kh  for  the  transient  solution  to  produce  damping  at  resonant 
frequencies. 


7.1  Illustration  of  time-dependent  trapped  wave  behavior 

The  diagram  below  (Figure  23)  illustrates  schematically  the  time-dependent  trapped  wave 
behavior  at  a  single  wave  number.  Given  time  t  (sec),  Mx-M2-  round  (t/2T )  where  the 
propagation  time  T  from  the  ground  to  the  turning  point  is:  T  -  mnL/kx  U0  (sec)  and  M, 
an  integer  function  of  time,  is  the  trapped  wave  round-trip  counter.  In  Figure  23,  za  is  the 
desired  analysis  level  and  the  “0”  subscripts  denote  surface  values.  The  hatched  region 
depicts  the  region  between  the  analysis  level  and  trapping  level  or  turning  point  at  z, ,  and 
tfr  is  the  time  elapsed  for  the  wave  in  the  present  trapped  wave  cycle. 


Fig.  23.  The  time-dependent  trapped  wave  domain  for  a  single  mode 


If  (za  <  zt)  and  M,  =  0 ,  the  trapped  wave  influence  is  not  included  and  fjtr  =  0.  Now, 
M{=M2  +  \  if  the  wave  front  lies  within  the  hatched  zone,  that  is,  if  the  wave  front  is 
crossing  za  on  its  way  up  to  the  turning  point  or  if  it  has  been  reflected  from  the  turning 
point  prior  to  reaching  za  on  its  way  down.  It  is  important  that  the  round-trip  counter  for 
the  trapped  waves,  M?,  be  incremented  after  each  round  trip.  Like  the  steady-state 
solution,  the  MWFM-3  transient  solution  only  accounts  for  a  single  turning  point  in  the 
vertical  column  for  each  wave  number  mode. 

7.2  Discussion  of  results  for  MWFM-3  transient  model  comparison: 

Calculations  displayed  in  Figures  24,  25,  and  26  are  compared  with  Broutman  et  al’s 
(2006)  example  results  for  their  Fourier-ray  solutions.  The  input  parameters,  like  those 
used  in  Figures  16  and  22,  are  L  =  2.5  km  for  the  linear  zonal  wind  shear  profile 
definition,  a  =  2.5  km  for  the  bell-shaped  hill  radius,  ho  =  100  m  for  the  hill  peak,  N  = 
0.01  s'1,  and  a  surface  wind  speed  of  Uo  =  10  m  s'1.  Calculations  were  performed  on  a 
1024  x  1024  horizontal  grid  with  approximately  1  km  grid  spacing.  Calculations  used  to 
construct  the  vertical  cross-sections  were  perfonned  at  1  km  vertical  intervals.  Broutman 
et  al  (2006)  also  obtained  results  similar  to  those  of  MWFM-3  using  a  numerical  model 
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Figure  24.  Comparison  of  vertical  velocities  calculated  at  z  =  2.5  km,  t  =  3  hours  for  MWFM-3 
transient  solution 
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Figure  25.  East-west  vertical  cross-section  of  vertical  velocities  along  hill  latitudinal  center  line  (y  =  0) 
for  MWFM-3  transient  solution  at  t  =  1  hour. 
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time  =  3  hours 
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akin  to  the  one  developed  by  Lipps  and  Hemler  (1982).  The  numerical  model  had  an 
absorbing  sponge  near  its  top  and  lateral  boundaries. 


At  z  =  2.5  km  for  the  horizontal  plan  view  (Figure  24)  at  t  =  3  hours,  the  calculated 
vertical  velocity  magnitudes  (~  +  0.21  /-  0.30  m/s)  agree  relatively  well  with  the  article’s 
Fourier  (+0.23,  -0.19  m/s),  and  numerical  model  results  (+0.25,  -0.21  m/s)  (which  are  not 
shown),  respectively.  The  time-dependent  behavior  pattern  agrees  qualitatively  with  the 
published  figures.  However,  the  calculated  trapped  mode  wavelength  (~  1 1  km)  is 
significantly  shorter  than  the  article  figures’  wavelength  (--15  km).  Figures  25  and  26 
depict  the  vertical  velocity  as  a  function  of  altitude  and  distance  downstream  from  the  hill 
location  for  latitude  y  =  0  at  two  different  times:  t  =  1  hour  and  t=  3  hours,  respectively. 
The  zonal  vertical  cross-section  results,  compared  with  the  published  values  (Broutman  et 
al,  2006),  also  exhibit  a  shorter  wavelength. 

Initially,  a  problem  was  encountered  with  the  vertical  structure  of  the  eigenfunctions. 
Unless  a  small  imaginary  component  (~  5  x  10'7  radian/m)  was  added  to  the  horizontal 
wave  number  for  the  x-direction  kx ,  the  eigenfunctions  did  not  damp  with  height  above 
1 0  km.  Furthermore,  when  the  imaginary  wave  number  component  was  included,  the 
peak  magnitudes  dropped  to  ~  50%  of  the  non-damped  amplitudes.  This  problem  was 
overcome  by  defining  the/;,  component  in  terms  of  the  vertical  integral  up  to  the  analysis 
level  and  allowing  the  phase  component  always  to  be  non-zero  (see  Eq.  23  and  the 
following  forms).  The  appropriate  procedure  for  handling  the  phase  functions  was  not 
clear  from  the  literature  and  was  discovered  after  some  trial  and  error. 

These  results,  as  was  the  case  for  the  uniform  steady-state  solutions,  appear  sensitive  to 
the  phase  functions  or  vertical  integrals.  The  phase  function  fa  is  not  defined  explicitly 
in  the  journal  article  (Broutman  et  al,  2006)  and  is  referred  to  alternately,  in  other  papers, 
as  the  phase  integral  for  vertically  propagating  waves  or  as  the  phase  for  a  ray  incident 
upon  the  caustic.  The  source  of  figure  discrepancies  may  be  due,  in  part,  to  lack  of 
knowledge  of  the  precise  definition  of  the  phase  function  fa. 

8.  APPLICATION  OF  MWFM-3 

The  Hawaii  12  December  2002  case  used  previously  to  compare  Smith’s  Three  Layer 
model  results  with  balloon-based  measurements  of  vertical  velocity  is  also  used  to 
demonstrate  the  behavior  of  the  MWFM-3  model  for  a  specific  application.  In  Figure  27, 
the  top  two  plots  illustrate  the  MWFM-3  steady-state,  fully  discrete,  hydrostatic  solutions 
at  altitudes  of  5  km  and  10  km  above  sea  level  using  the  padded  Hawaii  terrain  associated 
with  Figure  7a.  and  background  wind  and  temperature  fields  similar  to  the  basic  state 
listed  in  Table  3.  The  arrows  on  the  plots  depict  the  horizontal  wind  direction  at  altitude. 
The  plots  are  centered  at  the  same  location  as  depicted  in  Figure  7a.  The  resultant 
MWFM-3  model  vertical  velocities  range  from  +  6.3m/sto-4.3m/sat5  km.  The  major 
axis  of  the  vertical  velocity  fields  or  bands  is  roughly  perpendicular  to  the  wind  direction 
at  5  km.  At  10  km,  the  MWFM-3  vertical  velocities  are  noticeably  weaker  with  a  range 
of  +  1.6  m/s  to  -  1.8  m/s.  Also,  the  major  axis  of  some  vertical  velocity  bands  is 
approximately  parallel  to  the  flow  direction  at  10  km. 


z  =  5  km  MSL  (hydrostatic) 
wmax  =  +  6.3  m/s.  Aw  =  0.2  m/s 
wmin=  -  4.3  m/s  (+ w  contours  dark) 


z  =  10  km  MSL  (hydrostatic) 
wmax  =  +  1.6  m/s,  Aw  =  0.1  m/s 
wmin=  -  1.8  m/s  (+ w  contours  dark) 


Broutman’s  transient  solution  for  12  Dec  2002  100  UTC  (grid  line  spacing  =  200  km) 
(personal  communication,  2007) 


Figure  27.  MWFM-3  results  using  “real-world”  Hawaii  terrain  and  atmosphere  for  12  Dec  2002 
The  basic  state  for  the  upper  plots  is  taken  from  GFS  data  valid  at  12  Dec  0600  UTC  while  the  lower 
plots  use  wind  and  temperature  data  from  Hilo  rawinsonde  launch  at  12  Dec  1200  UTC. 


Dr.  Broutman  (personal  communication,  2007)  applied  the  transient  wave  MWFM-3 
model  for  altitudes  up  to  18  km  for  the  same  location  using  the  Hilo,  Hawaii  1200  UTC 
rawinsonde  for  the  background  velocity  and  temperature  profile.  Transient  wave 
solutions  for  the  5  and  10  km  levels,  corresponding  to  the  steady-state  solutions  described 
above,  are  depicted  in  the  lower  two  plots  of  figure  27.  The  spatial  distributions  of  the 
steady-state  MWFM-3  vertical  velocity  solutions  (upper  panel)  are  similar  to  the  transient 
wave  solutions  depicted  in  the  lower  panels  using  the  Hilo  rawinsonde.  Differences  in 
lower  boundary  terrain  and  background  state  specification  may  contribute  to  the 
differences  in  vertical  velocity  magnitudes  between  the  two  MWFM-3  versions. 


9.  RECOMMENDATIONS 

The  application  of  Fourier  transforms  to  practical  physical  problems  has  become  more 
routine  with  the  advent  of  modem  small  computers.  The  mountain  wave  models,  Smith 
Three-Layer,  and  Broutman  MWFM-3,  are  based  on  a  strong  theoretical  foundation 
which  has  been  developed  over  many  years.  This  foundation  is  well  documented  in  the 
literature.  Both  Smith  Three-Layer  and  Broutman  MWFM-3  have  undergone  limited 
verification  using  aircraft  and  satellite  data  (Smith  et  al,  2002  ;  Eckermann  et  al,  2006). 
These  models  are  worthy  of  testing  in  an  operational  prediction  environment,  while 
remaining  cognizant  of  their  physical  limitations. 

Recently,  the  Navy  has  modified  the  Smith  model  so  that  it  can  be  applied  to  multiple 
layers,  and  allows  both  critical  layers  and  turning  of  the  wind  using  arbitrary  two- 
dimensional  terrain  fields  for  the  lower  boundary.  The  modified  model  is  run  using 
operational  data  at  the  Naval  Research  Laboratory  (Doyle  and  Jiang,  2006)  for  regions 
such  as  the  Sierra  Nevadas  and  Rocky  Mountains  and  for  altitudes  up  to  12  km  above  sea 
level.  The  Smith  Three-Layer  model  has  not  been  employed  toward  prediction  of  waves 
in  the  stratosphere  above  1 5  km. 

The  Broutman  MWFM-3  model,  as  described,  has  been  developed  to  include 
tropospheric  and  stratospheric  wave  propagation.  It  has  no  theoretical  limitations  on  the 
number  of  model  levels  or  the  top  height.  The  model  results  in  the  literature  were 
typically  displayed  to  30  km  above  sea  level.  This  high  altitude  characteristic  favors 
MWFM-3  testing  in  an  operational  stratospheric  forecast  environment. 

The  Air  Force  Weather  Agency  (AFWA)  can  implement  MWFM-3  in  experimental  form 
in  regional  windows  using  operational  NWP  model  data  for  the  background  states. 
Implementation  details  can  be  provided  by  Dr.  Broutman  to  numerical  modelers  at 
AFWA.  Dr.  Broutman  can  advise  AFWA  modelers  on  how  best  to  apply  the  MWFM-3 
model.  MWFM-3  model  output,  based  on  operational  forecast  models,  can  be  posted  on 
the  Joint  Air  Force  &  Army  Weather  Information  Network  and  compared  with  other 
experimental  mountain  wave  forecast  models  including  MWFM-2  (Marks  and 
Eckennann,  1995)  and  the  S-Layer  Turbulence  Model  (Sinclair  and  Kuhn,  1991). 
MWFM-3  model  output  can  also  be  verified  using  stratospheric  in-flight  observations  as 
they  become  available. 
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APPENDIX  A.  EIGENVALUE  ANALYSIS  OF  THE  VERTICAL 
STRUCTURE  EQUATION 


A  Fourier-transformed  steady-state  elliptic  partial  differential  equation  for  vertical 
velocity  w  (m/s)  representing  the  vertical  structure  equation  for  a  two  dimensional  system 
(x,  z)  is: 

|^  +  (/2-C)w  =  0  (Al) 

where  w(k,z)  —  wr  +  iwt ,  (i  =  V=T)  is  the  complex  vertical  velocity  transform  or 


eigenvector: 


w(k,z)=  i — |  w(x,z^e'k'dx  (m  s"1  rad"1) 


\f2n  •*“' 

z  is  height  or  altitude  (m)  ,  k 2  is  the  square  of  the  one-dimensional  horizontal  wave 
number 


(A2) 


k  -  —  ( rad  / m)  (A3) 

Ac 

with  Lx  being  the  wavelength  in  the  direction  (x)  of  the  horizontal  wind  velocity  at  1  km 
height,  and 

f2/  \  N2  1  d2U  .  2.  , 

l~(z)  =  — r - —  ( rad  /nr)  (A4) 

y  U  U  dz 

is  the  Scorer  parameter.  N  is  the  Brunt-Vaisala  frequency  (s'1)  and  U  (m/s)  is  a  horizontal 
wind  speed.  The  hat  symbol  A  designates  those  variables  defined  in  the  transfonn  or 
wave  number  domain.  The  expression  for  the  vertical  dispersion  relation  or  wave 
number  is  written  as 

m2  -  I2  -k2  ( rad  !  m1 ) 


The  Scorer  parameter  / 2  is  calculated  using  a  vertical  profile  of  temperature,  pressure, 
relative  humidity,  and  the  horizontal  wind  speed  and  is  defined  as  a  real  function.  The 
eigenvalue  k2 is  treated  as  complex.  In  the  analysis,  the  eigenvector  amplitude  w  =  wj e' 

is  described  by  its  modulus  |vv|  and  phase  angle  (p  where 

|w|  =  ( w2  +  vv2  )A  (m2/radian  sec) 

<j>  =  tan'1  ( wt  /wr )  (radians) 


The  solution  to  (A  1)  as  a  boundary  value  problem  is  formulated  here  as  a  general 
complex  matrix  for  a  system  of  homogeneous  linear  equations: 

Ax-b  (A5) 

where  x  is  the  complex  eigenvector,  A  is  the  operator  matrix,  which  contains  both  the 
differential  operator  and  dispersion  relation,  and  b  is  the  right-hand  side  term.  Equation 
(A5)  is  solved  for  the  eigenvector  x  by  specifying  the  wavelength  or  eigenvalue  k2 
repeatedly  over  a  range  of  wavelengths.  Each  eigenvalue  is  associated  with  a  column 
eigenvector  with  dimension  equal  to  the  matrix  row  dimension  or  number  of  vertical 
levels.  The  wavelength  or  eigenvalue  which  produces  the  maximum  eigenvector 
modulus  is  selected  as  the  “dominant”  wavelength  or  eigenvalue  mode. 


For  this  boundary  value  problem,  the  boundary  conditions  are  defined  (following  Shutts, 
1997): 


w(k, 0)  =  1  at  the  lower  boundary  z  =  0 

dw(k,zt ) 


dz 


-  +  imw 


(k,zt)  =  0  for  the  upper  boundary  z  =  zt 


(A6) 

(A7) 


This  upper  boundary  condition  is  often  referred  to  as  the  radiation  boundary  condition. 


To  formulate  the  linear  operator  matrix  A,  Eq.  (Al)  is  approximated  by  the  finite 
difference  fonn: 


Wj+l-2Wj 


+  w, 


Az1 


■  +  m,jWj  =0 


(A8) 


where  A z  is  the  height  differential  separating  evenly-spaced  vertical  levels.  Equation 
(A8)  possesses  a  second-order  numerical  error  associated  with  truncation  of  the  Taylor 
series.  Similarly,  the  upper  boundary  condition,  Eq.  (A7)  is  approximated  in  forward 
finite  difference  form  as: 


w.  ~  wi- 1 

’  —  +  im/.wj=0 


A z 


(A9) 


which  has  a  first-order  numerical  error.  Given  the  lower  boundary  condition,  Eq.  (A6), 
the  matrix  solution  starts  effectively  with  level  j  =  2.  Grouping  terms,  level  2  is  written: 


w-,  -  2  w. 


V-,  ^  Wry  2  *  f\ 

- r— —  +  m2w2  =  0 

A z 


(A10) 


From  levels j  =  3  through  levels  j  =  n-  1 ,  where  n  is  the  top  level,  Eq.  (A8)  can  be 
expressed  as: 


wJ+l  +  w  ,  2  2 

— - yJ—+  m - T  \  w,  =0 

Az2  {  J  A z- 


(All) 


For  level  j  =  n,  using  the  upper  boundary  condition  (A9): 


im ,  H - 

1  Az 


w.  . 

w, - ^  =  0 

1  Az 


(A  12) 


The  operator  matrix  A  is  constructed  to  solve  for  the  column  eigenvector  f  from  levels  j 
=  2  through  levels  j  =  n.  This  yields  a  tridiagonal  matrix  of  order  (n  - 1)  x  (n  - 1) .  The 

problem  can  be  solved  numerically  for  x  using  a  subroutine  solver,  given  the  general 
tridiagonal  matrix  A  and  right-hand  side  column  vector  b,  or  using  a  relatively  simple 
non-iterative  algorithm  akin  to  Gaussian  elimination  for  solving  second  order  ordinary 
differential  equations  as  a  two  point  boundary  value  problem  (Lindzen  and  Kuo,  1969). 
Both  methods  yield  identical  results  for  the  dominant  mode.  For  the  calculations 
described  below,  the  dominant  eigenmode  was  selected  as  the  wavenumber  that  has  the 
maximum  eigenfunction  modulus  over  of  a  range  of  200  wavenumbers  ranging  in 
wavelength  from  0.25  km  to  50  km. 


An  advantage  of  the  eigenanalysis  relative  to  Fourier  wave  methods  appears  to  be  that  it 
allows  relatively  high-resolution  continuous  vertical  profiles  of  static  stability  and  wind 
speed  or  shear  (Scorer  parameter)  to  serve  as  the  basic  state  for  determining  a  vertical 
profile  solution  using  the  matrix  inversion  method  required  by  vertical  discretization. 

The  eigenanalysis  also  permits  quick  determination  of  the  dominant  wavelength  modes 
associated  with  the  basic  state  and  their  vertical  structure.  A  disadvantage  is  there  is  no 
provision  for  direct  inclusion  of  terrain  in  the  lower  boundary  condition  in  this  approach. 

A.l  Sample  calculation 

The  eigenvalue  calculations  are  sensitive  to  the  specification  of  the  Scorer  parameter, 
especially  to  whether  the  second  tenn  in  the  right  hand  side  of  Eq.  (A4)  —  the  wind 
curvature  term  —is  retained,  but  appear  less  sensitive  to  the  number  of  model  levels  or 
matrix  dimension.  The  Scorer  parameter  is  calculated  following  the  numerical 
procedures  described  by  Shutts  (1997)  using  a  height  differential  Az  =  100  meters  and  the 
horizontal  wind  speed  component  parallel  to  the  wind  direction  at  z  =  1  km.  This  altitude 
corresponds,  presumably,  with  the  mountain’s  ridge  top  altitude.  Test  calculations  (not 
shown)  were  performed  for  selected  horizontal  wavelengths  and  varying  numbers  of 
matrix  levels  using  a  sample  thermosonde  profile.  The  results  for  individual  modes 
exhibited  sensitivity  to  the  inclusion  of  the  second  tenn  in  the  Scorer  parameter  (the  shear 
derivative)  in  the  lower  10  km.  When  the  matrix  dimensions  were  increased  from  200  to 
300  levels,  little  change  was  observed  in  the  eigenfunction  behavior. 


Figure  At:  Eigenvector  profiles  for  the  dominant  mode  (#  72,  Lx  =  18  km)  and  the  Scorer  parameter  profile, 
without  the  second  term  on  the  right  hand  side  of  Eq.  (A4). 


The  features  of  the  dominant  mode  for  a  sample  rawinsonde  profile  (rokwn201),  at  a 
wavelength  of  18  km  (mode  72),  are  plotted  as  a  function  of  altitude  in  Figure  Al.  These 
include  modulus,  phase,  and  the  Scorer  parameter  with  the  second  term  involving  the 
second  derivative  omitted.  The  modulus  peaks  around  8  to  9  km  with  a  phase  angle  of 
180  degrees.  This  implies  a  zero  or  negligible  imaginary  component  for  the  eigenvector 
over  the  region  of  its  maximum  amplitude. 


Behaviors  of  the  dominant  modes  for  1 8  rawinsonde  profiles  as  a  function  of  wavelength 
are  illustrated  in  figure  A2.  The  peak  moduli  for  the  dominant  modes  were  detennined 
objectively  while  the  dominant  mode  tops  were  detennined  by  inspection  of  moduli 
contours  plotted  as  a  function  of  wave  number  and  altitude.  These  results  imply  the 
depth  of  the  dominant  mode  is  inversely  proportional  to  the  wavenumber  magnitude,  or 
proportional  to  the  wavelength.  Shallow  dominant  modes  are  associated  with  shorter 
wavelengths  while  the  taller  or  deeper  dominant  modes  are  associated  with  longer 
wavelengths.  In  addition,  the  logarithm  of  the  wavelength  of  the  dominant  modes  is 
somewhat  proportional  to  the  logarithm  of  the  eigenmode  moduli. 


Figure  A2:  Dominant  eigenvector  mode  tops  (km)  and  peak  moduli  (m2  s"1)  as  a  function  of  wavelength  (km). 
Each  open  circle  represents  the  dominant  eigenmode  for  a  single  Scorer  parameter  profde. 


References 


Lindzen,  R.S.,  and  H.L.  Kuo,  1969:  A  reliable  method  for  the  numerical  integration  of  a 
large  class  of  ordinary  and  partial  differential  equations.  Mon.  Wea.  Rev.,  97,  732-734. 

Shutts,  G.,  1997:  Operational  lee  wave  forecasting.  Meteorol.  Appl.,  4,23-35. 


APPENDIX  B.  REPRESENTATIONS  OF  THE  AIRY  FUNCTION 


Airy  function:  Ai(z)  =  —  £  cos(}t3  +  zt^dt 


(integral  representation) 


1.  Using  ascending  series  for  small  \z\ : 

Ai(z)  =  cj(z)-c2g(z)  ci  =  0.355028053887817  c3  =  0.258819403792807 

f(  .  /  x  y  3 V*+1  r(2/3  +  k) 

~  0  (3k)\  r(l/3)  ’  8^Z’~  V(3A:  +  1)!  r(2/3) 

k  =  1,  2, 3, .  .  .  r\  =  T  (r  +  1)  z  can  be  real  or  complex:  r  is  real 


2.  Using  asymptotic  expansion  for  large  \z\ : 


1  QO 

Ai(z )  — /r~^z_>Vf]r(-l)  where 

2  n 


r  2  } 

£  =  — z 
3 


r(3*  +  j) 
5Akk\T{k  +  \) 


d/(-z)  ■ 


-  n  '-z  • 


f  /zA 

sin 

l  A) 

3.  Auxiliary  functions  for  large  negative  real  z: 

Ji(-z)  =  z"^[/1(^')cos^'  +  /2(4')sin4']  (seep.  All  of  Abramowitz  and  Stegun,  1972) 

4.  Auxiliary  functions  for  large  positive  real  z  =  x  +  Oi: 


Ai(z)-^z'^e  cf{~  C)  (seep.  475  of  Abramowitz  and  Stegun,  1972) 


Figure  Bl.  Airy  function 
along  the  real  axis 


Fast  parameterization  of  Airy  function  where  z  is  real: 


1.  z  >1.0  :  use  auxiliary  function  for  large  positive  z  (item  4.) 

2.  -0.5  <  z  <  1.0  :  use  polynomial  fit 

3.  -7.0  <  z  <  -0.5  :  use  Lighthill  (1978)  approximation  for  substantial  negative  z 

Ai(z)  ~  n  ^  cos(f 

4.  z  < -7.0  use  auxiliary  function  for  large  negative  z  (item  3.) 
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Figure  B2.  Ai(z)  =  jAije'^  Airy  function  in  the  complex  plane  ( x ,  iy) 
z  —  x  +  iy  modulus  black  contour  lines 
phase  ( ° )  red  contour  lines 


ACM  machine  algorithm  (Amos,  1986)  Ai(x  +  iy)  =  \Ai\e“l 


Ascending  series  (e.g.  Miller,  1946)  Ai(x  +  iy )  =  \Ai\e"t’ 


Figure  B3.  Comparison  of  machine  and  series  solutions  for  the  Airy  function  in  the  complex  plane  (x, 
iy).  Contour  analysis  follows  the  same  format  as  figure  B2. 


